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Abstract 

We investigate the effect of excluded volume interactions on the electrolyte distribution around 
a charged macroion. First, we introduce a criterion for determining when hard-core effects should 
be taken into account beyond standard mean field Poisson-Boltzmann (PB) theory. Next, we 
demonstrate that several commonly proposed local density functional approaches for excluded 
volume interactions cannot be used for this purpose. Instead, we employ a non-local excess free 
energy by using a simple constant weight approach. We compare the ion distribution and osmotic 
pressure predicted by this theory with Monte Carlo simulations. They agree very well for weakly 
developed correlations and give the correct layering effect for stronger ones. In all investigated cases 
our simple weighted density theory yields more realistic results than the standard PB approach, 
whereas all local density theories do not improve on the PB density profiles but on the contrary 
deviate even more from the simulation results. 
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I. INTRODUCTION 



Understanding the behavior of charged macroions in solution is an importantproblem in 
fundamental science [jj as well as in industrial P and biological applications (J Charged 
stabilized colloidal dispersions are present in paints, inks, pharmaceutical products and are 
used in the fabrication of nanostructured materials 0, 5], Q] . These systems serve also as a 
primitive model for the crowded cellular environment that represents numerous biomacro- 
molecules and cellular polymers j?l gj- What all the applications above have in common is 
that when a charged macroion is immersed in an electrolyte solution, it is surrounded by 
counterions to balance the surface charge. The charged macroion surface along with the 
neutralizing diffuse layer of counterions is usually refereed as the electric double layer, the 
understanding of which is crucial for describing the behavior of such systems. For instance, 
the stability of colloidal dispersion depends on the distribution of small ions around the 
colloid. The electrophoretic mobility of the solution also can be rationalized in terms of 
the ion diction jmflQand naost of the electa! reactions oeeo r in this 
interfacial region 

As a result, there has been a considerable effort to describe the density profile around 
the macroion for different macroion geometries. The earliest theory that had significant 
success was the Poisson-Boltzmann (PB) approach. Its versions for planar geometry, the 
so called Gouy-Chapman theory jl^l l^ . can be solved exactly. It also has an analytical 
solution for an infinitely long linear macroion confined to a cylindrical cell Q, ^| , whereas 
only a numerical solution can be obtained in the case of a spherical geometry. The major 
flaw of this mean-field approach is that it neglects all correlations between the ions. For 
a long time, integral equation theories have been developed to adequately describe dense 
systems of electrolytes, and recently field theories have become very popular in calculating 
correlation corrections to the mean field PB approach, see i.e. Refs. jjlQQl for overviews. 
However, since the treatment of size effects is mixed with the electrostatic correlations, in 
many approaches it becomes difficult to identify the role of each effect. And finally integral 
equation theories work well at high densities when excluded volume contributions are very 
strong, whereas they are problematic in the low density regime. 

It would therefore be desirable to have a theoretical framework which retains the sim- 
plicity of the early attempts, but also accommodates correlation effects - something that 
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can be done within density functional theories. It is possible to rigorously rewrite the par- 
tition function of, say, a system of charged colloids, as a density functional jscj, in which 
the contribution beyond mean field is included as an additive correlation correction to the 
free energy density. The functional form of this correction is unknown and one has to use a 
reasonable Ansatz for it. The spirit is very similar to the fundamental problem of integral 
equations, where one also has to make an educated guess (namely, the closure relation). 
However, in the case of a functional this involves a free energy density expression rather 
than a relation between two- and three-point functions. It thus relies on a different kind of 
intuition and thus permits some complementary insight. 



A number of density functional prescri 
lations into account have been proposed (21 



;ions for taking both size and electrostatic corre- 



22, 



23, 



24j . These theories are able to reproduce 




to some extent the density profile of charged systems. However, since they treat both size and 
electrostatic correlations together, the origin of the result is not clear. Recently we adopted 
a different approach. We studied systems of point-like counterions (therefore no size effects) 
and addressed the question of when the electrostatic correlations become relevant. For treat- 
ing these correlations we proposed the Debye-Huckel-Hole-Cavity functional [25 . 
theory relies on a Debye-Hiickel treatment of the One Component Plasma (OCP) 
in which the short-distance failure of linearization is overcome by postulating a correlation 
hole. Since beyond a certain density the resulting OCP free energy density is a concave 
function of density, this favors the development of inhomogeneities. In the pure OCP these 
are balanced by the homogeneously charged background. However, if one uses the OCP 
free energy density as a correlation correction to the mean-field functional describing the 
double layer at a charged surface, one has all the charge opposite to the counterions located 
on that surface, rather than homogeneously distributed as a stabilizing background. The 
consequence is that the double layer becomes unstable and all ions collapse onto the surface, 
an effect which has been termed "structural catastrophe" 0, E^] • To prevent this effect we 
introduced a spherical exclusion region where no background can be found. The prescrip- 
tion for finding the size of such an exclusion serves both to keep the theory self-consistent 
and to establish the range of validity of the PB approach. Comparisons of the ionic charge 
distribution around a charged cylinder and a charged sphere showed a very good agreement 
with the simulations for both monovalent and trivalent counterions [26^ . 

Having studied how to take into account the electrostatic correlations, we address in this 
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paper the relevance of excluded volume correlations. We present a validity criterion for 
the PB approach by constructing a parameter whose value approximately indicates when 
excluded volume correlations are expected to become relevant. We then test several local 



density approaches, that have been advocated [32|, |33j , demonstrating that they all fail to 
take the size correlations into account, and that they even lead to an instability of the solution 
beyond a certain ion size. In order to circumvent this, a weighted density functional based 



on a simplified Tarazona approach 



34. 



35 



36[ is introduced. Our results are compared with 



Monte Carlo (MC) simulations, showing very good agreement for the cases of moderately 
developed hard core correlations, and even for strongly electrostatically interacting systems 
in both zero salt and non-zero salt cases. 

The remainder of the paper is organized as follows. In Sec. II we discuss validity of the 
PB approach for a colloidal system with non-point like counterions and show how the size 
effect can be incorporated into the model. Details of the used numerical methods are given 
in Sec. III. These are followed by the results and discussions presented in Sec. IV, and our 
conclusion in Sec. V. 



II. SIZE CORRELATIONS WITHIN DENSITY FUNCTIONAL THEORY 

Consider a spherical colloid of radius r c and negative charge Z, which is located in the 
center of a spherical cell of radius R c . This cell represents a bulk colloidal solution with 
colloid volume fraction = (r c /R c ) 3 . The counterions are taken as positively charged hard 
spheres of diameter a and valence v and N = Z/v of them provide the neutrality of the cell. 
The solvent is modeled as uniform dielectric background of dielectric constant e, and the 
strength of the electrostatic interactions is defined by Bjerrum length 

lB = A ^1, T ' W 

where q is the unit charge. In the non-zero salt case, N s positive and N s negative salt 
ions are also included. Here we assume that all ions have the same size and valence as the 
counterions. The average charge distribution is described by local densities n_(r) for the 
coions and n + (r) for the counterions and positive salt ions. These are defined for r Q < r < R, 
where tq = r c +a/2 is the distance of closest approach between the macroion and the particles 
and R = R c — a/2 (See Fig. Q). Therefore, the volume fraction of the bulk electrolyte </> e is 
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FIG. 1: Colloidal cell model. The system is defined by five independent parameters: number of 
counterions, N, number of salt ions, N s , and three characteristic sizes: tq, R and a measured in 
units of Ibv 2 . See Appendix for more details. 



given by 

a 3 (N + 2N.) 

*' = m-ri) ■ ( ) 

The effective surface charge density a should be defined in terms of r rather than r c , i.e. 
a = Z/(Anrl). This will be used later when defining the plasma parameter of the system. 

The central task of a density-functional theory is to derive an analytical expression for 
the free Helmholtz energy functional that upon minimization gives the density profiles of 
the free ions in solution. Its simplest form is given by the Poisson-Boltzmann functional, 
namely 

F PB = J d 3 r |A; B Tn + (r)[m(n + (r)A 3 ) - l] + k B Tn_(r) [in (n_(r)A 3 ) - l] + f e i(r)^) 

which includes the translational entropy of the ions (A is the thermal de Broglie wavelength) 
and all electrostatic interactions represented by / e ;(r). The electrostatic interactions within 
the mean-field approximation are given by 

/ezO) = vq[n+(r) -n_(r)]^(r) , (4) 



where ip{r) is the total electrostatic potential created at position (r) by the fixed macroion 
and all ions together. The minimization of Eq. (jSJ) with respect to n + (r) and n^(r) gives 
the Boltzmann density distributions 

n± (r) = n° ± e*^ (r) , (5) 

where parameters n° + and n°_ are defined by the charge neutrality condition. In spherical 
geometry, Eq. |5] together with the Poisson equation 

V 2 ip{r) = — - — [n+(r) — n-(r)\ (6) 

and the boundary conditions at r = ro and r = R comprises a fully defined Poisson- 
Boltzmann problem. The problem with this approach is that the ions are considered as 
point charges in some average electric field and both electrostatic and excluded volume 
correlations between them are not taken into account. When do these correlations matter? 
Comparing the PB predictions to simulation results, one has found out, that, for point-like 

evant when the plasma parameter, T 2 d = A/vrcr/^f 3 , 



ions, electrostatic correlations become re 
becomes larger than 1 



25 



26 



37 



38, 



39l . Here we address the question under what 
condition the excluded volume effects become significant. For a uniform hard sphere liquid 
we know that the radial distribution function changes from monotonically decaying to non- 
monotonic at volume fractions of the order of ~ 0.2, which will be the reference volume 
fraction in our further estimates. In the case of confined liquids, ions tend to concentrate 
close to the surface and their concentration at the surface, especially in charged systems, 
may be much higher than that in the bulk j^lj]. Therefore, to access the hard-core effects 
one should consider the average volume fraction within the first layer of counterions close 
to the colloid surface. Since the colloid is much larger than the ions it can be approximated 
as (for simplicity the salt is not included) 

2 ra 

4> s = —7- / n PB (x)dx , (7) 

D JO 

where npn(x)£ 3 = T^/ti^xT^ + l) 2 is the exact solution obtained in the planar geome- 



try 



42|. Here x = x/£, where I = Ibv 2 . After integration we find that the volume 



fraction close to the macroion is given by 



o 3 r 4 

(8) 



3(2rl d a + i) 



where a = a/£. Then, for <p s < 0.2, there are only weakly developed excluded volume 
correlations, and the PB approximation should be still valid. For larger values of <p s we 
expect to see some layering effects close to surfaces. This criterion is strictly valid only 
for a planar geometry, but is expected to approximately hold for sufficiently large spherical 
or cylindrical macroions where the curvature effects (oc 1/R) are negligible. An analogous 
formula which takes the curvature into account can also be derived for a cylindrical PB cell 
for which the contact density is known analytically. Since electrostatic correlations were not 
taken into account, we do not expect this simple analysis to hold for T 2 d ^ 2. Beyond this 
value, the force- distance curves between charged plates cease to be monotonic, and beyond 
T 2 d — 2.45 attractions even between like charged macroions can occur [37]. These effects 
are results of correlations between different double layers (like, for instance, ion interlocking 
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43j) that are stronger than the size effects we are describing here. Note that, Eq. (jHj) is 
designed to give a limit of validity of the mean-field approach. For high ionic radius <p s can 
become larger than 1, loosing its relation to the actual volume fraction in the system. 

Correlations can be included into the PB model by adding to the PB free energy Fps an 
excess free energy term, 

F = F PB + F ex . (9) 

The excess free energy F ex originates from internal interactions within the system and it is 
unknown. In principle, it can account for both the hard-core repulsion and the electrostatic 
correlations. Since we want to test the excluded volume effects within the range of ionic 
strength in which they overcome the electrostatic correlations, the latter will be neglected 
within our approach. 

There are a number of various functionals which can be used to include hard-core effects 
in uniform liquids of a given density n. Within the local density approximation (LDA) 
one can choose one of these free energy density expressions, f e x[n], and replace the uniform 
density by a local one so that the total excess free energy reads 

F ex = k B T J d 3 rn(r)f ex [n(r)} . (10) 

For our system we can take n(r) = n + (r) + n_(r) which means that the hard-core effects 
are treated identically for both positive and negative ions. The idea behind Eq. (IK)jl is quite 
simple. A particle at r is supposed to be affected by only the particles around it, in a range 
given by the interaction. If the range of the interparticle interaction is much smaller than 



the typical length for variations in n{r), the system can be divided in small subvolumes of 
nearly constant density and each of them can be treated as part of a homogeneous system. 
If we take ions as charged hard spheres, we can use the free energy density derived, for 
example, from the Carnahan-Starling (CS) equation of state [44], namely 

0(r)(4-30(r)) 



fcsW)] 



(11) 



(1-0W) 2 

where <f>(r) = it a 3 n(r)/6 is the volume fraction occupied by the free ions. For denser liquids, 
the accuracy might be improve d by using the more precise virial expansion for the Percus- 
Yevick theory for hard spheres |45{], namely 



fvir[(j>(r)) = |40(r) + 50(r)^ + 6.120(r)' j + 7.O20(r) 4 + 7.9O5(/.(r) 5 + 9.42O80(r) b | .(12) 

The last two expressions agree up to second order in local volume fraction 4>(r). 

A simple form of F ex can also be derived from the free volume (fv) expression for a lattice 
gas (3^, |3| , namely 

" 1 



fvl 



k R T / d 6 r 



— n r 



In (l - n(r)a 3 )) 



(13) 



where a is the lattice spacing. With this form of functional the excluded volume effects can 
be explicitly incorporated into the PB equation. However, its density expansion is different 
from that of hare sphere. Another expression based on the free volume concept which can 

'(r) 



be found in Ref. 



fv2 



-k B T I d a rn(r)\n I 1 - 



(14) 



gives lower order density terms similar to Eq. (jllj) and Eq. (|12p. 

The assumption of smooth variations of n(r) being within the characteristic range of 
hard-core interactions is only valid for sufficiently small ionic radii a. Consequently, as we 
will show later, all the functionals above underestimate the densities close to the colloid and 
completely fail to give the correct density profile in the more interesting cases where stronger 
variations in n(r) are observed. Moreover, all expressions above have a singularity at certain 
value of volume fraction. This reflects the fact that the bulk density cannot exceed some 
upper limit. If a local density is higher than that, for example due to a charged surface, the 
local density functional does not converge and no density profile can be obtained. 

In order to circumvent this difficulty a number of weighted density approaches (WDAs) 



have been developed 



H 35,y, 



47, 



48, 



for the description of neutral hard sphere 



solutions, and some rather complex methods have already been proposed for charged sus- 



pensions 



21 



22, 



The prescription we have followed here, in the spirit of the generalized 
van der Waals theory of Nordholm and co-workers 50] , is to introduce the non- locality of the 
free-energy density functional through a coarse-grained density distribution n(r), which at 
each point is a non-local functional of the local density n(r). This can be pictured as a mean 
density around point r averaged over a volume related to the range of the interactions. In 
this context, the local density in Eq. (fTU|) is replaced by some weighted density n(r), namely 

F ex = k B T J d 3 rn(r)f ex [n(r)} (15) 



where 



nir 



J d 3 r'w(\r - r'|)n(|r'|) . (16) 



The weight function w(\r — r'|) should be chosen to give reasonable direct correlation func- 
tions which are functional derivatives of F ex [n]. The most important are the first- and 
second-order correlation functions, defined as 

- W (i7) 

(2)/ n S 2 F ex [n] 
& >{r,r ) - 



8n(r)8n(r') 

In Tarazona's [3, H| approach one assumes that the weight itself is also density 
dependent and can be expanded in powers of the weighted density as follows 

w(|r — r'|) = wo(r) + wi(r)n(r) + W2(r)n(r) 2 + ... . (18) 

If we substitute this expression into the direct correlation function in Eq. ([17)1 . the resulting 
expansion can be set equal to the direct correlation function of a uniform hard sphere 
fluid [5^. This way it is possible to obtain the weight function that up to second order is 
given by Q Q, Q 
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Ana 3 L 

Wl (r) = [0.475 - 0.648^ + 0.113 (-) ], r < a (20) 

= [0.288- -0.924 + 0.764-- 0.187 (-)\ a < r < 2a 
r a \aJ 

= 0, r > 2a 

, , 57ra 3 r r /r\ 2 n . 

«fc(r) = lii-[ 6 - 12 - + 5 (J ] ' r<a (21) 
= r > a. 

Since our aim is not to precisely describe the hard-sphere effects but just to access their 
relevance, we will employ the simplest form of the weight function that is the first term in 
Eq. El or a constant weight given by Eq. We will refer to this weight as WD AO 
whereas the weight defined by Eqs. [T9"ll2*Tlwill be called WDA2 and will be used to validate 
our results. For a pure hard sphere fluid, WDA0 reproduces the discontinuity in the direct 
correlation function predicted by Percus and Yevick However, it overestimates 

the range of the correlation function, especially at high densities, when com par ed to the 



n tunction, especially at mgn densities, wnen compare 
HQ or to direction dependent weights 22, 



density dependent weig hts [2U 124 

or to direction dependent weig hts 22, 49, 51J. Having; 
this in mind, we will concentrate on the systems for which size plays a relevant role but the 
differences between the constant weight and more sophisticated approaches do not affect our 
main conclusions. 

Measuring all lengths in units oil = Ibv 2 reveals that the full partition function of our cell 
model depends on five system parameters. The observables, for example, reduced density 
profile n(r) = n(r)£ 3 , remain constant under rescaling which does not change the following 
quantities: the number of counterions N = Z/v, the number of salt ions N s , the reduced 
distance of closest approach f = r /£, the r /R ratio, and the reduced ion diameter a = a/£ 
(see Appendix). The same holds for PB theory and is also true for our WDA theory. For 
the WDA it implies that both the WDA free energy correction and the weight function have 
to obey this restriction. 
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III. NUMERICAL METHODS 



In this Section we give details of the numerical methods used to study the cell system 
described in Sec. II. Three different ways were employed to find the ion distribution in the 
cell. The first one was a direct Monte Carlo simulation of the cell model which gave us 
reference data to test the theoretical results. The density profile minimizing a given free 
energy functional was obtained using numerical iteration until it converged to the equilibrium 
charge distribution. Another way of minimizing the functional was by Monte Carlo sampling. 
Some technical details of these three methods are summarized below. 

A. Monte Carlo simulation. 

Within this approach we simulate the cell model exactly as it is - all ions are taken as 
charged hard spheres of diameter a confined between two spherical shells of radii r c and 
R c . To gather the statistics of charge distribution the ions are moved around the cell and a 
single ion move is either accepted or rejected according to the usual Metropolis probability: 



where AE is the difference between the system energy after and before the move. Since 
the density profile is highly anisotropic, a combination of two types of moves was found to 
improve the efficiency of the sampling. An ion was either inserted at a random position in the 
cell or randomly displaced within a cube centered at its current position. The former allowed 
for efficient exploration of low density regions, whereas the latter proved to be efficient close 
to the colloid where a successful insertion of an ion could be a rare event due to the high 
packing fraction. The frequency of using one or the other move as well as the displacement 
range were adjusted to give an about 50% acceptance rate. 

B. Iterative functional minimization. 

When minimizing a functional containing a non-zero correlation term, Eq. El becomes 
dependent on the excess chemical potential 




(22) 




(23) 
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and reads 

n± (r) = n «e T ^^ (r) -*± w . (24) 

Once the expressions for n<± are derived for a given functional, we can find the ion distribu- 
tion which satisfies both Poisson equation and Eq. |^ Integrating the Poisson equation over 
a spherical shell of radius r, and using the Gauss theorem, an integro-differential equation 
for the electric field E(r) can be obtained. Consequently, the optimum density profile can 
be obtained from the numerical iteration of this equation until convergence is achieved. 



C. Monte Carlo functional minimization. 



Within this approach the ion position is described only by its distance from the colloid, 
r, and this uniquely defines the density distribution, n(r). Each MC step consists of moving 
an ion to a new trial position r$ < r < R. This move is either accepted or rejected with 
probability Q 

7i = min{l, exp (-0AF)} , (25) 

where AF is the free energy difference after and before the move and it is explicitly given 
by the functional we minimize. This method was found to be more stable and worked much 
faster than the iterative procedure, though the final result did not depend on the numerical 
approach used. 



IV. RESULTS AND DISCUSSION 



In this Section we compare how well the different density functional approaches described 
in Sec. II capture the excluded volume interactions. First, we apply numerical techniques 
described in Sec. Ill to two colloidal systems already studied in the literature j&J. Then 
we perform a systematic analysis of hard-core effects by investigating a number of different 
systems with and without added salt. 



We start from considering two salt free systems which were also used in |33( to study the 
excluded volume effect in colloidal solutions. In both systems r$ = 50A, R = 100A, a = 10A, 
and lb = 7 A. The number of monovalent ions is different and it is N = 200 in system (a) and 
N = 500 in system (b). Note that already for 200 ions some packing effects are expected to 
be seen because (p s = 0.25. Figure [21 shows the ion density distribution close to the colloid 
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(r-r )/a (r-r )/a 

FIG. 2: Ion distribution close to the colloid surface measured in systems with r$ = 50 A, R = 100A, 
a = 10A, Ib = 7 A, and containing (a) N = 200 and (b) N = 500 monovalent ions. Each curve 
corresponds to a particular method used: MC is the result of the Monte Carlo simulation; Fes, 
F v ir, Ffvi and Fj v2 are obtained using LDA with a form of the excess free energy given by Eas-ITO 
I141 correspondingly; WDA0 is the constant weight curve and PB is the numerical solution of PB 
equation without any hard-core corrections. 



obtained using different approaches for both (a) and (b). Comparison between the PB and 
the MC curves for system (a) reveals that hard-core repulsion decreases condensation by 
pushing ions away from the colloid. This effect is captured well by the WDA, whereas it is 
significantly overestimated by all LDAs - the predicted contact densities are too low. The 
highest packing fraction achieved in LDA calculations for system (a) was below the critical 
value and the equations converged. 

In the case of 500 ions, the LDA with functional given by Eq. was found to be nu- 
merically unstable. The iterative functional minimization failed to converge, while some 
convergence could be still achieved by explicitly limiting the highest density at a -3 in the 
MC sampling of the functional. The result, however, depended on the number of bins and 
therefore was not physical. A plateau close to the colloid, as seen in j^, was observed 
only for relatively large bin sizes, while smaller bins resulted in a saw-like density profile 
(not shown here). The other LDAs, while still converging, overestimated hard-core effects 
and also resulted in unphysical profiles. Moreover, all local approaches missed the layering 
clearly captured by WDA0 and observed in simulations at distance a from the colloidal 
surface (Fig. [2Kb)). 
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TABLE I: Parameter cf) s of different ionic sizes (columns), Bjerrum (rows) lengths and plasma 
parameter (rows). Both Bjerrum lengths Ib and ion sizes a are given in units of ro- Packing effects 
are expected to be seen in those systems for which <j) s > 0.2. Unrealistically high values of 4> s 
observed for large ions indicate inapplicability of PB and failure of LDA for these systems. 



Below we consider a model system in which parameters r and R = 5r are kept fixed, 
whereas the ion size and the Bjerrum length are varied. The cell containing N = 100 
monovalent ions of size ranging between 0.1r < a < 0.8r is studied at four Bjerrum 
lengths 0.1r , 0.2r , 0.3r and 0.4r . These correspond to plasma parameter 0.5 < T 2 d < 2.0 
which enables us to investigate hard-core effects in systems with weak, moderate and strong 
electrostatic correlations. The parameter <p s calculated for each of these 32 systems and 
given in Table predicts that, at all Bjerrum lengths, the size correlations are expected to 
be seen for ion sizes a > 0.3ro, where (f) s is already greater than 0.2. The Monte Carlo 
data show that, indeed, the ion density profile is concave for a = 0.1r and a = 0.2ro but 
develops a convex region at distance about a from the colloid surface for larger ion sizes at 
all four Bjerrum lengths. This indicates some packing taking place which is well captured 
by our s -criterion. Figure EH shows both the reference Monte Carlo and WDA0 density 
profiles calculated for different ion sizes at fixed Bjerrum length l B = 0.2r (T 2 d = 1-0). 
The development of layering with increasing the ion size is well captured by the non-local 
functional approach, whereas none of the LDAs exhibits any layering and therefore are not 
shown in Fig. El 

Another way of checking how well correlations are captured by a particular excess free 
energy functional is to compute the osmotic pressure. In real systems this pressure also 
depends on correlations between ions of different cells, something which is not taken into 
account within the cell model approximation. So by pressure we refer to the pressure exerted 
on the rigid wall at r = R of our cell model. Within the simulations, the pressure is given j^j 
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FIG. 3: Ion distribution close to the colloid surface measured in systems with plasma parameter 
T2d = 1-0. Odd (a) and even (b) ionic diameters are shown separately for sake of clarity. Only 
WDAO curves are marked while all MC curves are presented as solid lines and can be identified by 
the corresponding closest dashed WDAO curve. All distances are measured in ionic diameters so 
the formation of the second layer of counterions is always expected to be around 1. 

by the contact density at r = R: 

IT = k B Tn(R). (26) 

For the density functional approach, this exact expression should be corrected to recover 
the free energy functional after integrating pressure over the volume. The correction term 
is generally small and for simplicity we will directly compare contact densities predicted 
by different methods. Figure 0] shows both the colloid contact density n(ro) and boundary 
density n(R) given by the Monte Carlo simulations and different local and non-local density 
approaches for systems at fixed Bjerrum length Is = 0.1r (T 2 d = 0.5) as a function of the 
ionic diameter a. The colloid contact density is informative of how well a certain method 
works at the most packed region of the system and can also be related to the pressure. 
Figure Efb) shows that the local approaches underestimate n(ro) at high ionic radius when 
compared to the simulations. We found that even for small ionic sizes PB density profile was 
closer to the MC reference data than the results of any of the locally "improved" functionals. 
Moreover, the limitation of the local approach is illustrated by the failure of the LDA's to 
converge at large a (no data are shown for large ions). Both WDAs we used here give 
consistent results when compared to MC simulations. 

Now, we concentrate on the case of relatively large ions of diameters a > OAro for which 
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FIG. 4: Contact density (a) n(ro) and boundary density (b) n(R) predicted for systems with 
^2d = 0-5 by different methods as a function of ionic diameter. All theories at a = are identical 
to PB which slightly underestimates the contact density due to ignoring electrostatic correlations. 
This effect reduces as the ion size increases. The legend is organized similar to that of Fig. |2] with 
one more line WDA2 corresponding to the weight by Tarazona given by Eqs. I19B21I There is a 
slight positive difference in n(R) predicted by WDA2 and WDAO not seen clearly in (b) at this 
scale. 



layering is clearly observed. According to the s -criterion, the density profiles of such systems 
should deviate from PB and the hard-core interactions can be even more significant than 
the electrostatics in some cases. Figure El shows the integrated ion fraction 

i r 

P(r) = — / dr 4vrr 2 n(r) (27) 



N 

for systems with ionic diameters fixed at (a) a = 0.4ro and (b) a = 0.6ro for plasma 
parameter T 2 d 0.5, 1.0 and 2.0. Clearly, a larger plasma parameter leads to an increased 
condensation (the curves are shifted up) - an effect which is governed by the electrostatics 
and also present in PB theory. One could expect that for high T 2 d electrostatic correlations 



would be a dominant effect 26] , and significant deviations to PB theory and also to our WDA 
corrected DFT should arise due to electrostatic correlations, which we did not account for. 
However, under the investigated circumstances, the ionic size plays a more relevant role. The 
hard core effects lead to packing effects that overcompensate the electrostatic correlations. 
For a = 0.4ro, the density profiles are well captured by WDAO, which is always much closer 
to the MC data than to the PB result (not shown here). However, for a = 0.6r , the structure 
of packing becomes important at T 2 d = 2.0. For this system, the WDA2 weight improves the 
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FIG. 5: Integrated fraction of counterions obtained using Monte Carlo (MC) and one constant 
weighted density approach (WDAO) for plasma parameter I^d = 0.5, 1.0, 2.0 and for ionic diameters 
(a) a = 0.4ro and (b) a = 0.6ro. In (b), the WDA2 curve is also shown for the case of r 2 d = 2.0. 



result, showing that the deviation between the simulations and the constant weight WDAO 
is not due to the electrostatic correlations but rather to hard-core effects. Beyond this point 
a more sophisticated functional should be used to capture the local packing. 

In principle, the addition of salt can lead to new correlations due to ion-ion correlations 
and screening. For high electrostatic salt couplings, T s = |, where d is the distance of closest 
approach of ion and coion, ion clusters can also appear, that change the ion distribution 
considerably However, this effect can be overcome by the hard core if the ions are 

large enough, rendering T s < 1. Below we consider the systems from Tab. |J with added 
salt. Two cases of N s = 10 (10% of salt) and N s = N (100% of salt) are studied to 
represent moderate and high amount of salt. Since addition of the salt ions into the cell 
would increase the packing fraction, here, we prefer to keep it constant adjusting accordingly 
the cell radius R. This is partially justified by the fact that the ion distribution close to 
the colloid weakly depends on the cell size for our system parameters. Figure |B1 shows both 
the positive and negative charge density profiles obtained using simulations, WDAO and 
PB for the system with T 2 d = 1.0, a = 0.4ro with 10% of salt. Due to the screening, the 
effect of electrostatic correlations is less profound than in the zero-salt case. The agreement 
between the simulations and WDA is therefore improved at higher salt concentrations and 
lower plasma parameter. The integrated charge profiles for several systems employing our 
highest plasma parameter T 2 d = 2.0, are shown in Fig. The WDA captures the same 
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FIG. 6: Density profiles of positive and negative ions obtained using MC simulations, WDAO and 
PB for the system with T 2 d = 1-0, a = OAtq and 10% of salt. 




FIG. 7: Integrated fraction of counterions obtained using (a) Monte Carlo (MC) and (b) one 
constant weighted density approach (WDAO) for different ion sizes at plasma parameter T2d = 2.0 
and different amount of salt (N s = 0, 10% and 100% of N). The sizes vary from top to bottom as 
a = O.lro, a = 0.4ro and a = 0.6ro, 



trend in the charge distribution as provided by the MC data. The presence of salt does not 
change the layering as is observed for this ionic radius when the salt is not present. One 
should be aware that things will be more complicated if one considers asymmetric salt (in 



valence and size), or 



can occur 
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arge ion sizes, since then more complicated effects like overcharging 



18 



V. CONCLUSION 



In this paper we studied the effects of adding various local and non-local free energy 
functionals to the PB free energy functional to include the effects of an ionic hard core. We 
started from calculating a layer volume fraction S using the PB approximation, that gives 
a criterion for recognizing when size effects become relevant and the simple PB approach 
has to be modified. We tested this criterion for a system consisting of a charged spherical 
colloid and its counterions confined to a spherical cell, and studied a number of parameter 
combinations where the PB approximation fails. For including size correlations, four local 
and two non-local density functional approximations were employed. The local theories 
were always found to overestimate the hard-core effects, creating an exclusion region close 
to the colloid for large ionic radii. Beyond a certain ionic radius, all the considered LDAs 
diverge and produce meaningless results. The failure of the LDA is also captured by the 
increasing divergence between the LDA and the MC contact densities which is seen when 
the ionic radius is increased, and the absence of any layering effect in the LDA. Due to this 
observations we note that the inclusion of the LDA correction into PB actually worsens 
the agreement of PB with simulation results. On the other hand, we demonstrated that a 
simple weighted density approach for the excluded volume interaction was able to capture 
the main features of the ionic density profile. The introductions of a more sophisticated 
weighted density approximations such as Tarazona approach jjyj improves the agreement 
with the simulation, but it does not bring any new physics to the problem. If some salt is 
included, under certain parameters the main effect is the increase of screening of the electro- 
static correlations. Therefore, the system can be adequately described by the PB approach 
supplemented with an excluded volume WDA. More complicated effects are expected to 
appear at sufficiently high plasma parameters and higher salt concentrations. To treat those 
within density functional theory, a combination of hard-core and electrostatic correlations 
along the lines of Refs. (3, [y] will probably be required. 
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APPENDIX 



The canonical partition function Z of the colloid surrounded by its counterion in a cell 
model is given by: 



N+N a N s 



/i v TJV S i v s 
nn 

1=1 7 = 1 



d 3 pid 3 rid 3 pjd 3 rj 



hW+Xs)h 3N s(N + N s )\N s 



(28) 



where N = Z/v is the total number of counterions, 2 N s is the total number of positive and 
negative ions of salt. The Hamiltonian H = T + V splits into kinetic and potential degrees 
of freedom. In the classical description employed here the kinetic part T will contribute the 
usual factor \- ZN - &N s to the partition function, where A is the thermal de Broglie wavelength. 
The potential energy can be expressed as 
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where the first two terms are related to the electrostatic interactions and the last is re- 
sponsible for the hard-core repulsion. The specific form of this term is not relevant here. 
After rescaling all length by £, i.e. introducing f := r/£, the total partition function can be 
rewritten as 
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where d = a/£. 
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In this form it becomes evident that appropriately scaled thermal observables like the 
integrated charge density (measured in units of £~ 3 ) or the pressure (measured in units of 
/csT£ -3 ) are invariant under system changes which keep the number of counterions N, the 
number of salt particles N s , the rescaled colloid size f = r /£, the rescaled ion radius a, 
and the volume fraction <fi constant. 

Poisson-Boltzmann theory shows the same invariance property, as does the approximate 
density functional theory we are proposing in this paper. 
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